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Glass transition where viscosity of liquids increases dramatically upon 
decrease of temperature without any major change in structural 
properties, remains one of the most challenging problems in con- 
densed matter physics [1, 2] in spite of tremendous research efforts 
in last decades. On the other hand disordered freezing of spins in a 
magnetic materials with decreasing temperature, the so-called spin 
glass transition, is relatively better understood [3, 4]. Previously 
found similarity between some spin glass models with the structural 
glasses [5, 6, 7, 8, 9] inspired development of theories of structural 
glasses [10, 11, 12, 13, 14] based on the scenario of spin glass tran- 
sition. This scenario though looks very appealing is still far from 
being well established. One of the main differences between stan- 
dard spin systems to molecular systems is the absence of quenched 
disorder and the presence of translational invariance: it often as- 
sumed that this difference is not relevant, but this conjecture is 
still far from being established. The quantities, which are well de- 
fined and characterized for spin models, are not easily calculable for 
molecular glasses due to the lack of quenched disorder which breaks 
the translational invariance in the system and the characterization 
of the similarity between the spin and the structural glass transition 
remained an elusive subject still now. In this study we introduced 
a model structural glass with built in quenched disorder which alle- 
viates this main difference between the spin and molecular glasses 
thereby helping us to compare these two systems: the possibility of 
producing a good thermalization at rather low temperatures is one 
of the advantages of this model. 



Glass transition 
metry Breaking 



Random Pinning | Static Correlation Length | Replica Sym- 



Dramatic slowing down of relaxation process in almost all 
liquids when supercooled below the melting temperature 
still lacks a proper explanation [1, 2]. The increase of relax- 
ation time in deep supercooled regime is so impressive that 
it becomes extremely difficult for the system to reach equi- 
librium in experimental time scales and eventually the liquid 
falls out of equilibrium with further decrease of temperature 
to under go a calorimetric glass transition. This transition 
is defined at a temperature where the viscosity of the liquids 
becomes 10 13 poise [18, 19]. It is clear that this transition is 
ad hoc in nature and depends crucially on the choice of the 
parameter, but the main question, which remains to be an- 
swered, is whether there is a true thermodynamic transition 
below this calorimetric transition temperature. 

Many approaches to understand this remarkable slowing 
down in the dynamics of supercooled liquids invoke the ex- 
istence of a cooperative length scale [20] associated with the 
collective rearrangements of particles and its divergence at 
the elusive glass transition. This scenario of glass transition 
is very similar in spirit to the critical phenomena seen in the 
continuous phase transition. The slowing down is believed to 
be caused by the difficulty to rearrange sets of ever increas- 
ing number of particles in a collective fashion with decreasing 
temperature or increasing density. Attempts to estimate this 
cooperative length scale remained one of the major difficulties 
due to the lack of identification of an order parameter char- 
acterizing the structural glass transition. This is primarily 
due to lack of growth of an obvious order in the system with 
decreasing temperature. 

Progresses made in recent years to identify such a length 
scale is really encouraging. Dynamic heterogeneity length 
scale from the analysis of 4-point density-density correlation 



function [22, 21, 23, 24], Point-to-set length scale [25], patch 
length scale [26], length scales associated with non affme dis- 
placement of particles [27], from finite size scaling of config- 
urational entropy [24] and density of states [28] are a few to 
name. Unfortunately there is still no general consensus about 
the importance of these length scales to glass transition and 
their relation to each other. The main hurdle in reaching 
such a goal is that these length scales are only accessible in 
computer simulation studies because of the requirements of 
microscopic details to compute them. So these length scales 
can only be estimated in small parameter range where it grows 
very modestly thereby making it almost impossible to see the 
divergence while approaching the glass transition. 

Based on these ideas of growing length scale and re- 
markable similarity seen in the dynamics of p-spin glassy 
model with the structural glasses, Kirkpatrick, Thirumalai 
and Wolynes [10, 13] proposed the Random First Order the- 
ory (RFOT) of glass transition. This theory is very much 
in spirit with the Adam-Gibbs Theory [16] proposed much 
earlier. This theory seems to suggest that p-spin glass mod- 
els and structural glasses should belong the same universality 
class. However there is an inherent difference between these 
two models, that is the existence of quenched disorder in spin 
glass and not in the structural glass models. Lack of quenched 
disorder in structural glasses makes it rather difficult to cal- 
culate the quantities like spin-glass type order parameter and 
susceptibility: moreover the whole low temperature phase can 
not be accessed by simulations. In a recent mean field and 
renormalization group study [29], it was proposed that with 
random pinning one can explore the ideal glass phase as in 
the temperature - concentration of frozen particles plane there 
exists a critical point where relaxation time does not diverge 
while going from liquid to glass phase. Similar studies done in 
the Mode Coupling Theory frame work [30, 31] also confirm 
this picture. All these studies and some other recent studies 
[32, 33, 28, 34, 15] clearly show that exploring the glassy state 
in the random pinning geometry can be very insightful with 
the added advantage of built in quenched disorder. This can 
enable us to do a replica theoretic calculations for these kind 
of particle models to shed more light on the soundness of the 
replica approach to structural glasses. 

In general a system of particles of equal size interacting 
via a radially isotropic pairwise potential will have liquid to 
crystal phase transition with decreasing temperature for di- 
mensions d < 3. At higher dimensions it was shown that 
crystallization is strongly suppressed [35]. So to form a glass 
it is very important at least for dimension smaller than 4 to 
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Fig. 1. Top Left Panel: Pair correlation function g(r) for TV = 512 system size for different temperatures. The curves for different temperatures are shifted vertically 
for clarity. Top Right Panel : Self part of the two point correlation function Qs(t) for different temperatures. The inset shows that a-relaxation time T a calculated as the 
time when Q s (t) goes to 1/e of its initial value. The line is the fit to the data using Vogel-Fulcher-Tamman Formula (VFT) with Tk ~ 0.17. The dynamic transition 
temperature (MCT) is estimated to be ~ 0.394 ( red dashed line ). Bottom Left Panel : The extracted static length scale from the replica overlap correlation functions 
( see the method section for details ). The red circles, green squares and blue diamonds are from the system with N = 512, ./V = 2000 and TV = 4000 respectively. 
The estimation of the length scale from ./V = 512 system size seems to have some finite size effects which almost disappear for TV = 2000 and TV = 4000 system sizes. 
The black line is the fit to the data using the form £(T) ~ \T — T^p, with chosen Tx — 0.17 and ■y ~ 0.50 and the red line with dot is the same fitting where we 
allowed all the parameters to vary and the resulting T/ v ' 0.38 and 'y ~ 0.25. Inset shows the dependence of the relaxation time with the length scale for ./V = 4000 
system sizes. One can see some degree of universality in the relation between relaxation time and length scale (see text for details). Bottom Right Panel: Time-Temperature 
superposition for the Q s (t) for all system sizes N = 512, 2000, 4000 for all the studied temperatures. The very nice collapse of the data confirms that these model indeed 
has all the usual features of a glassy system. Inset shows the temperature dependence of T a for different system sizes. 



introduce frustration in the system to prevent it from quickly 
falling in to the crystalline global minimum. In spin glasses 
the random interactions between the spins are the source of 
this frustration and in structural glasses compositional disor- 
der, for example different sizes of the particles or the asym- 
metric interaction between different types of particles, usu- 
ally generates the required frustration. Due to presence of 
these extra degrees of freedom in general glass models they 
are often very hard to equilibrate at lower temperature as one 
needs also to equilibrate these extra degrees of freedom. It 
would be nice to have a glassy model system without extra 
degrees of freedom. In this article we propose to generate 
the required frustrations by random pinning: we have study 
a system where we do not have any randomness in the inter- 
action potential nor do we have any compositional disorder. 
Our model system consists of particles of equal size with some 
fraction pt mp of them frozen randomly in space. If sufficient 
number of these particles is frozen randomly in space at some 
high temperature there is enough frustrations in the system 
to force the system to remain in disorder state. As there 
is no other degrees of freedom to equilibrate apart from the 
positions of the particles we expect this model can be equili- 



brated to lower temperatures than the usual models to study 
the lower temperature phase of the supercooled liquid. With 
the quenched disorder, this model is also very attractive to 
compare it with the spin glass models. At high density of 
quenched particles no crystallization is present, however the 
glass phase is also absent: we will show that there is an in- 
termediate region with a small, but not too small fraction of 
quenched particles, where a part of the usual structural glass 
phenomenology survives and the system does not crystallize, 
also at very low temperatures. Our construction differs from 
that of [29] as far as the quenched disorder is crucial to avoid 
crystallization and finding a phase transition. 
Results : It turns out that one needs to freeze around 
pimp = 9% of particles to get a system which will not show 
crystallization at least for bigger system sizes (A'" > 250). 
For smaller systems this amount of frozen particles is found 
to be not sufficient to prevent the crystallization. With 
pimp — 11%, even the smaller system sizes do not show any 
tendency to crystallize. So for the present study we choose 
to work with pi mp = 11% for all the system sizes studied. 
We restricted ourselves to small system sizes mainly because 
we wanted to achieve full equilibration of the system to very 
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Fig. 2. Left Panel: The collapse of distribution of energy P(E) for different temperature on the distribution of P(E, To) with To = 0.425 in the Parallel tempering 
run for ./V = 100 system size according to the Eq.12 to check the equilibration. The nice collapse confirms that the very good equilibration is achieved using the Parallel 
tempering method for lower temperatures. Right Panel: The temperature dependence of the energy for different system sizes. One can see that the finite size effect is not very 
strong here. Inset shows the specific heat calculated from the fluctuation of potential energy for different system size. Some finite size effect can be seen here. 



low temperatures. The studied system sizes are in the range 
TV 6 {100, 4000}. In the top panel of Fig. 1, we showed the pair 
correlation function for different temperatures to confirm that 
there is no sign of incipient crystallization in this temperature 
range. In top right panel we showed the self part of the average 
two point density correlation function (Q s (t)) (see Dynamics 
in the Method Section for details) for TV = 512 system size 
for different temperatures in the range T £ {0.400, 1.000}. 
One can see the nice development of the plateau in the cor- 
relation function with decreasing temperature. In the inset 
of this panel, we have shown the temperature dependence of 
the a-relaxation time and the line is the fit to the data us- 
ing Vogel-Fulcher-Tamman (VFT) formula with the estimated 
divergence temperature Tk ~ 0.17 for TV = 512 system size. 
The dynamic transition temperature or the mode coupling 
cross over temperature is estimated to be Ta ~ 0.394 by a 
power fit. One should keep in mind that these extrapolated 
estimation of the divergence temperatures may not be very 
reliable as the range of the data is not very big. 

As our model is not translationally invariant due to the 
presence of the quenched disorder defining overlap between 
two replicas becomes easy and unambiguous. For example 
if we have a system with translational symmetry then while 
defining the overlap between two replicas we need to take in to 
account the fact that two replicas can be similar even though 
they may be translated or rotated in space [39]. We defined 
the overlap between two replicas using the window function 
w(x) which is 1 if x < 0.30 and zero otherwise. We say the 
two replicas are similar if for each particle at position r in 
replica 1 there is a particle of replica 2 within a sphere of ra- 
dius 0.30 (see Replica overlap in Methods Section for details). 
We also defined local overlap and calculated the corresponding 
spatial correlation function to extract the static length scale 
over which the replicas are similar. In the bottom left panel 
of Fig.l the extracted static correlation length of two replicas 
(see the Replica Overlap and Extraction of static length scale 
in Method section for details [37]). The finite size effect seems 
to die out quickly once we go to TV — 2000 system size for lower 
temperatures. The modest growth of the static length scale in 
this temperature range is very similar to other glass models. 
The inset shows the dependence of the relaxation time with 
this length scale (see the relation between length scale and 



relaxation time in Method section for details). The apparent 
universality for two different frozen particles density is also 
in agreement with the recent findings for usual glass models 
[36, 38]. 

To achieve equilibration at further lower temperature we 
implemented Parallel Tempering simulation methods follow- 
ing [40], but restricted ourselves to systems up to TV = 250. 
The details of the simulation method and the parameters used 
are given in Parallel Tempering methods parts in Method Sec- 
tions. We parallelized the Parallel Tempering method using 
Message Passing Interface (MPI) routines to speed up the sim- 
ulation. We run different replica in different computer cores 
and we found that the system can be equilibrated to very low 
temperature within reasonable CPU time (around 12 hours 
for TV = 250 particles using 16 replicas in 16 cores). We 
checked that system does not crystallize even at the lowest 
temperature studied. We believe parallel tempering methods 
works so well for our model as it has only positions to equi- 
librate and does not have any other compositional disorder 
which needs further equilibration. In left panel of Fig. 2, we 
have collapsed the probability distribution of potential energy 
P(E, T) for different temperatures according to the ansatz 
Eq.12 on the probability distribution of the reference temper- 
ature To = 0.425, to ascertain whether proper equilibration 
is achieved in our simulations (see Parallel tempering part in 
Method section for details). The nice collapse of the data 
indicates that the equilibration is achieved to a very good 
accuracy. The right panel of Fig.2 shows the temperature de- 
pendence of the average potential energy and the inset shows 
the corresponding specific heat calculated from the fluctua- 
tion of potential energy. Although average potential energy 
does not show strong finite size effects one can see somewhat 
strong finite size effect in the specific heat. The nicely de- 
veloped peak in the specific heat seems to be a precursor to 
a possible second order phase transition as also seen in [17]. 
The usual discontinuity in the specific heat seems to remain 
rounded in a volume independent faction. Notice that at the 
lowest temperatures we have simulated the specific heat drops 
to the DulongPetit value, suggesting that harmonic degrees of 
freedom are the mostly relevant ones at these temperatures. 

In Fig. 3, the temperature evolution of the distribution 
of overlap q is shown for six different system sizes. One can 
clearly see the change of the distribution from Gaussian to Bi- 
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Fig. 3. Temperature evolution of the Probability distribution of overlap q for six different system sizes N = 100, 150, 200, 250, 512, 4000. 
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Fig. 4. Left Panel : Temperature dependence of overlap for different system sizes. Bigger system size seems to show sharper change of overlap q with temperature. Middle 
Panel: The spin glass susceptibility is plotted as a function of temperature for the studied system sizes. The strong variation with temperature is really remarkable. Notice 
the strong system size dependence of susceptibility at lower temperature. Currently we can not say whether this susceptibility will diverge with decreasing temperature in 
thermodynamic limit as the system sizes studied here are still quite small. One needs to do bigger system sizes to find out the possible divergence of the spin glass susceptibility 
at lower temperature. Right Panel: Scaling collapse of the spin glass susceptibility using the scaling ansatz Eq.4. The exponents are rj ~ —0.1, V — 1.0 with as shown 
in the inset. 



modal with decreasing temperature. The distribution seems 
to deviate from the Gaussian one at temperature close to the 
temperature where the specific heat also shows peak as a func- 
tion of temperature for that system size. The shape of the dis- 
tribution as function of the temperature recalls what happens 
in mean field theory in the replica approach: below the crit- 
ical temperature a peak at higher values of q appears while 
the low q peak has an intensity that is proportional to the 
temperature. 

In left panel of Fig. 4, the temperature dependence of 
the overlap ( see Method section for detail definition ) for 
different system sizes are shown and one can clearly see 
that for larger system sizes the overlap seems to change 
more sharply. The corresponding susceptibility is given by 
Xsg — (N — N imp )[(5q 2 }], where (.) is the thermal averaging 



and [.] means the averaging over the different realizations of 
the disorder. Ni mp is the number of impurity particles in the 
system which is given by Ni mp = pi mp N. In middle panel of 
Fig.4, the susceptibility shows dramatic variation with tem- 
perature and at lower temperature the susceptibility seems 
to change quite strongly with system size also. For system 
size iV = 250 the susceptibility changes almost by a factor 
of 50 compare to its high temperature value. With this data 
divergence of the susceptibility at lower temperature in the 
thermodynamic limit can not be established but the strong 
increase of its value with the system size is encouraging, we 
have tried to collapse the susceptibility data for different sys- 
tem sizes using the finite size scaling ansatz (L = TV 1 / 3 ): 

Xsg = L 2 -^(l l '»\T-T k \), [1] 
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with rather poor results. A better collapse is shown in the 
right panel of Fig. 4 with r\ ~ —0.1, v ~ 1.0 and Tk{L) now 
being a function of system size L. It is important to notice 
that as Tk grows with system size it is reasonable that the 
divergence temperature will not go to zero in the thermody- 
namic limit. The data seems to indicate a possible non-zero 
temperature thermodynamical transition of second order in 
nature. The value of v is somewhat larger than the exponent 
v = 2/d argued in RFOT theory [10, 13] by scaling argu- 
ments. These results are very similar to the one obtained for 
finite range p-spin glass model in 3 dimension in [17]. The 
inability to collapse the data with a size independent critical 
temperature may be the effect of the strong finite size scaling 
corrections: the system sizes studied here are relatively small. 
A different scenario would be (in the thermodynamic limit) 
a strong increase of the susceptibility approaching the mode 
coupling transition temperature followed by a cross over to a 
different behavior at low temperature. We noted that the sus- 
ceptibility data for N = 4000, T >= 0.45 (that should have 
small finite size effects) are well fitted by a power law, with a 
critical temperature around 0.4, i.e. the putative mode cou- 
pling transition. Much more extensive simulations are needed 
to better understand the nature of the transition, if any. 

Mode Coupling Theory calculation in similar geometry for 
the hard sphere system in [31] suggests that at some critical 
density (pimp ~ 0.15 in [31]) of the frozen particles the time 
dependence of the two point correlation function will change 
from being two steps to one step. So to check whether we 
are not very close to this critical density we performed similar 
studies for different pi mp for N = 100 and calculated the over- 
lap distribution P(q). One can clearly see in Fig. 5, that there 
is no qualitative change in the behaviour with different pin- 
ning density. In right panel we show how < q > changes with 
Pimp as a function of temperature and the top right corner 
inset show the corresponding susceptibility. 
Conclusions: 

To conclude we showed how a very simple particle model with 
random pinning can be used to explore the glassy phase at 
deep supercooled regime and also showed how this model can 
be used to compute different spin glass correlation functions 
to shed light on the relation between spin glass transition and 
structural glass transition which is the basis for most of the 
recent theories of structural glasses. This model can also be 
used to find out the relation between different length scales 



as in this model the length scale is calculated directly from 
the spin glass order parameter correlation function. It would 
be extremely interesting to arrive to a precise determination 
of the phase structure of the model and to compare it with 
accurate numerical simulations. 



Materials and Methods 



Simulation details : 



4> {nj 



The interaction potential for Tij 

1 



< X c is given by 



[2] 



while it is for Tij > X c , where Tij is the distance between particle I and J, and 
X c is the dimensionless length for which the potential will vanish continuously up to 
q derivatives. The coefficients C2£ are given in [28]. 

Dynamics: We studied the dynamics using two point overlap correlation function, 
Q{t), defined as 

r N 

Q(t) = / drp(r,t )p(r,t + t ) ~ ^ w(\fi(t ) -fj(to + t)\) [3] 

i,J=l 

where p(f,t) is particle density at space point rat time t, and w(r) = 1, if 
T < d and zero otherwise. Average over the time origin to is assumed. The use of 
the window function [a = 0.30] is to remove the fluctuations in particle positions 
due to small amplitude vibrational motion. In Eq.3 the contribution only due to the 
self-term is denoted as Q s (t). The structural relaxation time T a is the time where 
Qs{r a ) = 1/e. 

Replica overlap and Extraction of static length scale: In this section we will 
explain how we have extracted the growing static length scale in this system. We 
followed the method used in [22]. The method is explained briefly below, we start 
with the following definition of overlap between two replica 



- N N 



N 



[4] 



i— 1 fc — i 

with W\X) = 1 for X < 0.3 else 0,. Now we define the coarse grain variable as 

fJ>(x) = y^S(x- Xi)qi, and defined f(r) = (fJ,(x) fJ,(y)) [5] 

i 

where T = \x — y\. Now to remove the natural oscillation in the function f(v) 
we divide the function by pair correlation function gir^j to define another function 
c ( r ) = f( r )/9( r )' and fitted the function to the following fitting function 

c(x) = a + exp(— x/£) [b + ccos(xd + e)] [6] 
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to extract the length scale We extracted the correlation length scale for the 
TV = 512, 2000 and 4000 system sizes which have been equilibrated using stan- 
dard molecular dynamics simulations using Berendsen thermostat [41] and we studied 
the system in the temperature window T £ {0.400, 1.000}. We averaged 
the data over 20 different realizations of the disorder. For last three temperatures 
T = 0.400, 0.420 and 0.450, we averaged the data over 40 realizations of the 
disorder. 

Relation between length scale and relaxation time : We start with the ansatz 
for the relation between relaxation time T a and the static correlation length 



r a (T) oc exp 



T 



[7] 



where Ao is a non-universal coefficient that depends on the details of the glass for- 
mer. Now at reference temperature T = To (the highest temperature in this case) 
we define the typical length to be Then the relaxation time at that temperature 



r a (T ) oc exp 
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So we have the following relation 
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[8] 



[9] 



As the pre-factor A(j is not known a-priori for different models we choose Ao = 
1.0, 1.11 for pimp = 0.090 and 0.110 respectively in Fig.l. We also choose 
\jj — 1.0 for these two cases. 

Parallel tempering methods : To equilibrate the system still at lower tem- 
perature, we have implemented parallel tempering method. We briefly mention 
the method here as details can be found in [40]. We construct a system con- 
sisting of TVf non-interacting subsystems (replicas), each consists of TV particles 
with phase space coordinate {T >l ,Q*}, where P % = {pi,P2 ■ ■ -piv} and 

= {Ql,Q2 ? ■ ■ ■ Qn} for the subsystem. The Hamiltonian of the i* 
subsystem is given by 



H i (P i ,Q i )=K(P i )+A i E(Q i ), 



[10] 



where A is the kinetic energy, E is the potential energy of the system and 
Ai G {-^1? • • • Am} is the parameter to rescale the potential energy. Now 
we perform molecular dynamics simulation of the whole system with Hamiltonian 
l~i = ■ Hi at a reduced temperature T = l//?0 using a modified isokinetic 
simulation method. Time step for the MD is taken to be St = 0.005. Now at each 
time interval of Atx = 1000(5i we exchange the parameter A between different 
replica I and J keeping all other things unchanged. The exchange is accepted using 
a Metropolis scheme, with a probability 



Wij = min(l,exp(— Aij)) 



[11] 



where Aij = f3o(Aj - A i )(E(Q t } - E{Q j )). We perform these steps for 
sufficiently long time such that we get proper equilibration of the system. This way we 
generate canonical distribution at inverse temperatures p% — \%pQ. Here we have 
also parallelized the code using MPI to speed up the simulation process. We have used 
Nl = 12 replicas between temperature 0.300 to 0.600 and another 12 replicas 
for temperature range 0.600 to 1.000 for TV = 100 and 150 system size. For 
TV = 200, we have used 16 replicas for 0.300 to 0.600 temperature range and 
12 between 0.600 to 1.000. For TV = 250, we only tried to equilibrate the 
system up to T = 0.330 and used 16 replicas between temperature range 0.330 
to 0.600 and 12 between 0.600 to 1.000. We have averaged the data over 400 
different realizations of the disorder for systems TV = 100, 150, 200 and 250. 
So the estimated computer time for TV = 250 system size is close to 10^ hours. 
To check the equilibration of the system we used the method as in [40] to rescale the 
Canonical distribution function using the formula 



PiiEiTj =A,/3o) = 



Pi(E)exp[(\j- X^PqE] 
fdE'Pi{E')exp[(Xi - Xj)f3 E'] 



[12] 



In equilibrium the left hand side of the above equation should be independent of % to 
within the accuracy of the data as can be seen in left panel of Fig. 2. 
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I. TIME-TEMPERATURE SUPERPOSITION FOR Q(t) 
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FIG. 1: Data for N = 512, N = 2000 and N = 4000. Left Panel: System size dependence of the relaxation function Q(t). One 
can see that there is very small finite size effect in these system sizes. Right Panel: The time-temperature superposition of Q(t) 
using the relaxation time r a obtained from the self part of Q(t) from the condition Q s {r a ) = 1/e. It clearly shows that the 
initial part of the relaxation function is very much same as that of the self part but deviates from it as the infinite time value 
of Q(t) is not in this case. The black solid line is the master curve one would obtain if one tries to do the time-temperature 
superposition for Q s (t) as done in Fig.[l] in main article. This is shown just to point out that the master curve for Q(t) also 
follows the same curve but deviates from it due to non-zero asymptotic value. Inset shows the temperature dependence of 
the infinite time value of Q(t) that is Q(oo) for all the system sizes. One can cleary see that there is hardly any system size 
dependence in this quantity. 
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In SMFigdJ we have shown the time-temperature superposition for the full overlap correlation function Q(t) with 
the same relaxation time r a which are used to collapse the data for self part of the correlation function, i.e. Q s (t). The 
figure clearly shows that initial decay of the correlation function is same as that of the self part only. Q(t) deviates 
from the master curve because of the non zero asymptotic value Q(oo). In this figure we also plotted the data for 
all system sizes. Different color indicates the temperatures and symbols distinguishes the different system sizes. The 
inset shows the temperature dependence of the asymptotic value of Q(t) i.e. Q(oo) and one can clearly see that there 
is very small finite size effects in this quantity. 

II. EXTRACTION OF STATIC LENGTH SCALE FROM THE REPLICA CORRELATION FUNCTION : 



As discussed in the method section of the main paper about how we extracted the length scale scale from the replica 
correlation function, we have presented here some data for the N = 4000 system size in SMFigJ^l We used the method 
explained in the Replica Overlap and Extraction of static length scale section of the Method section of the main paper. 
The nice fitting of the data over the whole range and for all temperatures indicates that the extracted length scale is 
very reliable. This can also be seen from the data shown in Fig.[l] of the main paper where we compared the results 
of this fitting for different system sizes. The data for N = 2000 is very close to the one obtained from N = 4000 
system size. This implies that this method of extraction of the length scale is very robust once we rule out the finite 
size effect in the data. 
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FIG. 2: Extraction of the static length scale for N = 4000 system size. We used the functional form c(x) = a + 
exp(— x/£) [b + ccos(xd + e)] to fit the data. The nice fitting over the whole range of the data indicates that the length 
scale extracted this way is very reliable. Data for different curve is shifted vertically by 0.01 from each other for clarity. 



III. CHECK OF EQUILIBRATION IN PARALLEL TEMPERING RUNS: 
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FIG. 3: Evolution of < q 2 > as a function of MD steps to check whether the time evolution of < q 2 > reached a plateau. We 
take the sample to be equilibrated if < q 2 > reaches a plateau value within the error bar. 



We performed two separate tests to check the equilibration of our sample : one is on the distribution of energy 
P(E) where we tried to rescale all the P(E) for different temperatures to collapse on the distribution of our reference 
temperature To as shown in the Fig. [2] of the main article using the ansatz Eq.[12] in Parallel tempering part of the 
Method section. Here we present the data for the second test. The test consists of checking when square of the replica 
overlap < q 2 > reaches a steady plateau value with time. In SMFigl3l we showed data for N — 100 system size and 
the data clearly shows that our simulation time is longer than the time required to reach a plateau value for < q 2 > . 
As we have started from two completely uncorrelated replicas their initial overlap is very small and over the time it 
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reaches the equilibrium value from below. 




FIG. 4: Comparison of susceptibility x f° r two sets of runs for different system sizes. Run2 is 2 times longer than Runf . Error 
bars are calculated by bootstrap method. 

In SMFigHl we have shown the spin glass susceptibility x calculated for all the 4 different system sizes for two 
different run length. Runl is 1 x I0 8 MD steps long and the Run2 is 2 times longer than the Runl. The value of the 
susceptibility is within error bar of each other. This also confirms that our data is well equilibrated. 



IV. SAMPLE TO SAMPLE VARIATION OF P(q) : 

In this section we have shown the sample to sample fluctuations of the distribution of P(q) for two different system 
sizes at the lowest temperature simulated for two system sizes N — 150 and 250. The SMFig. [5] shows P(q) for 8 
randomly chosen different samples. One can clearly see the strong sample to sample fluctuations of this distribution 
which is also very similar to spin glass models Some samples have only two peaks, with a deep valley in between, 
while other samples have a third peak. 



V. SKEWNESS AND KURTOSIS: 



In Fig. [6l we have shown the skewness (5) and the kurtosis (B) defined as 

s= mk> B = 3 -°-wSw with w 
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The usual Binder cumulant is equal to g(T) = B/2.0 according to our definition. Notice that g(T) is non-monotonic 
function of temperature. In conventional phase transition this is very useful quantity which has been used extensively 
in finite size scaling analysis for determining the critical temperature in the thermodynamic limit. This quantity for 
different system sizes cross each other at critical temperature because of the finite size scaling. 

At zero magnetic field in the Sherrington-Kirkpatrick model for spin glass and the short range Edward- Anderson 
model g(T) is always positive function of T and increases regularly with decreasing temperature and the skewness is 
identically zero. Unfortunately, as soon the magnetic field is different from zero already in the Sherrington-Kirkpatrick 
model (where mean field theory is correct) the behavior at the transition point for the cumulant is much more complex 
and also for 1024 spins we are very far from the asymptotic limit Q. 

The behaviour we see here is distinctly different from this simple case. Moreover it is non-monotonic and mostly a 
negative function of temperature over the almost whole range. This stems out from the fact that the distribution of 
q, P(q) starts to become non-Gaussian at quite high temperature. This behaviour is similar to the one seen for the 
finite range 3-spin model in 4 dimensions [3], where the skewness is strongly decreasing towards the large negative 
values when the temperature decreases. This similarity once again seems to suggest that this model may also in the 
same universality class that of short range 3-spin glasses which apparently have their own universality class. We need 
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further study to clearly find out the exact nature of the transition if any. 
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